Nucleation of colloids and macromolecules in a finite volume 
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A recently formulated description of homogeneous nucleation for Brownian particles 
in the over-damped limit based on fluctuating hydrodynamics is used to determine 
the nucleation pathway, characterized as the most likely path (MLP), for the nu- 
cleation of a dense-concentration droplet of globular protein from a dilute solution 
in a small, finite container. The calculations are performed by directly discretiz- 
ing the equations for the MLP and it is found that they confirm previous results 
obtained for infinite systems: the process of homogeneous nucleation begins with 
a long-wavelength, spatially-extended concentration fluctuation that it condenses to 
form the pre-critical cluster. This is followed by a classical growth processes. The 
calculations show that the post-critical growth involves the formation of a depletion 
zone around the cluster whereas no such depletion is observed in the pre-critical 
cluster. The approach therefore captures dynamical effects not found in classical 
Density Functional Theory studies while consistently describing the formation of the 
pre-critical cluster. 
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I. INTRODUCTION 



Recently, the process of diffusion-limited nucleation has been formulated within the con- 
text of fluctuating hydrodynamics^-. This approach gives a unified description of the forma- 
tion and growth of an unstable, pre-critical nucleus by means of thermal fluctuations as well 
as of the post-critical growth driven by a free energy gradient. In this dynamical approach 
to the description of nucleation, the fundamental quantity is a hydrodynamic field, the lo- 
cal concentration^. It can be used to derive the elements of Classical Nucleation Theory 
(CNT)2£ when the concentration-field is parameterized using the usual capillary approxi- 
mation that is implicit in CNT. More complex parameterizations are possible giving thereby 
a more fine-grained description of the nucleation pathway. By this means, it was found that 
the dynamical theory predicts an unexpected, non-classical pathway wherein nucleation is 
initiated by a long-wavelength, small-amplitude concentration fluctuation. Spatial localiza- 
tion of the mass in the fluctuation leads to the formation of a relatively large pre-critical 
cluster. This is followed by a second stage consisting of growth of the pre-critical cluster 
in accordance with more classical approaches (CNT, see, e.g., Ref. [f], Density functional 
theory (DFT), see, e.g., One advantage of the dynamical approach is that it 

can be used to determine the relative likelihood of different proposed nucleation pathways 
and it was found that the non-classical path is much more likely than those of more classical 
approaches, such as DFT, which inevitably show cluster formation as a spatially-localized 
processes^. 

All of the results described above were obtained for an infinite system. This leaves open 
the question as to whether the non-classical path will dominate in a finite system more 
typical of a real experiment. Since the physics of the non-classical path is rooted in mass 
conservation, it seems possible that in a finite system with finite mass, a different result 
could emerge. This could be particularly relevant in very small systems such as those typical 
of nano-fluidics, biological processes and lab-on-a-chip applications. A second question left 
open in previous studies is whether the use of parameterized profiles, no matter how complex, 
might be prejudicing the results. The purpose of this paper is to address both of these 
issues. Calculations based on the dynamical description of nucleation are presented for the 
nucleation of a bubble of dense protein-rich solution from a mother-phase of dilute protein 
solution in a small, finite volume. The calculations are performed by directly integrating a 



discretized version of the field equations that determine the most likely path in concentration- 
field space connecting the initial, nearly uniform dilute phase and a final state consisting of 
a large, stable droplet of dense solution within a small spherical cavity. It is found that the 
two-step mechanism is preserved in modified form and that further details of the process 
that were inaccessible with the parameterized models, such as the formation of a depletion 
zone during post-critical growth, are evident. Further, the additional resolution afforded 
by the direct discretization of the equations governing the nucleation pathway, as opposed 
to the use of parameterizations, allows to further distinguish the results of the dynamical 
theory from those of more standard DFT-based approaches. 

The next section presents the basic theory describing nucleation in terms of a temporally 
evolving concentration field. The stochastic description is discretized and the determination 
of the initial state, critical cluster, final state and most likely pathway between them is 
formulated. The third section presents calculations for the nucleation of protein-rich model 
globular protein from a dilute mother phase. The paper concludes with a discussion of the 
results compared to those obtained by other approaches and with some remarks concerning 
open questions. 

II. THEORY 

A. Stochastic evolution of the concentration 

The fundamental quantity in the present formulation is the local number concentration 
p(r). In the over-damped limit, the evolution of the local concentration is governed by the 
stochastic differential equation^ (SDE) 



where the noise, £ (r; t), is delta-correlated in both space and time. Even though this SDE 
involves a state-dependent noise amplitude, it can be showrt^ that it is Ito-Stratonovich 
equivalent so that there is no ambiguity in its interpretation. The functional F [p] is a 
coarse-grained free energy. In general, it is not equivalent to the free energy functional of 
DFT as the latter is an equilibrium quantity. However, it is generally assumed that the 
coarse-grained quantity is well-approximated by a mean-field approximation. Just to make 
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clear the connection, one expects that averaging Eq.flT]) over the noise would give something 
like the usual DDFT form 

a( " (r;( »-Sv ( pM)>vf^l (2) 



dt \r v , // . 5f) ^ 



p(r)^(p(r;t)> 



where, now, F [p] is indeed the equilibrium free energy functional of DFT. For spherically 
symmetric problems, the SDE can be reduced to 

dm(r;t) _ n ,_ 2 ^.^ d (5F[p\ 

p(r)->p(r;t) 

where the spherically symmetric noise, £ (r; t), is delta-correlated in both variables^. Here, 



D4nr 2 p(r;t)^-( 5 -fj^) + y/8nr*Dp(r;t)£(r;t) (3) 

dr \5p(r) J o(r) _> p(r . t) 



the quantity m (r; t) is the total mass within a ball of radius r, 

m (r; t) = Anr' 2 p (r'; t) dr'. (4) 
Jo 

For the free energy functional, a simple squared-gradient model will be used 

F[p] = J {/(p(r)) + ^(Vp(r)) 2 }dr (5) 

where / (p) is the bulk Helmholtz free energy per unit volume. Finally, it is interesting to 
note that the equation of motion can be written in the particularly simple and illuminating 
form 



dm (r; t) dm (r; t) 5F [p] 



™^i(r;t) (6) 



dt dr 8m (r) 

This shows that, in terms of the mass variable, the dynamics is gradient-driven with a 
fluctuation-dissipation relation. 



B. Discretization 
1 . Variables 

In general, SDE's of the form given here are only well-defined once they are discretized. 
Furthermore, discretization is necessary for practical calculations. Here, space will be dis- 
cretized via points = iA and the boundary of the finite volume of the container under 
consideration will be at r N = (N + 1/2) A. The masses at those points are m (r*) = for 
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< % < N + 1 and, of course, mo = 0. The concentration will be represented at intermediate 
points as 

P U + = p i+1/2 , 0<z<N (7) 
so that the relation with the mass is 

i-l 

m i = A 2^ 47rr j+l/2Pj+l/2 < > 47rr i+l/2Pi+l/2 = ^ ( 8 ) 

3=0 

Note that the exact linear relation between the masses and densities means that these 
variables are completely equivalent. Similarly, one has that 

d 1 / 1 d Id 



(9) 



drrii A4vr \rt_ 1/2 dp^ 1/2 r' +l/2 dp i+1/2 

With this discretization, the variable Pn+i/2 ls concentration at the wall and m^ + i is the 
total mass within the container. Since mass is conserved, the value of rriN+i is, by definition, 
constant and that of Pn+i/2 * s determined by mass conservation and the other densities via 
Eq.©. 



2. The free energy 

The free energy is discretized using the same integration scheme as relates the densities 
and masses, 

0F = A 4vrr| +1/2 /3/ ( Pj+1/2 ) + A ^ An-K r£ (10) 

Note that the upper limits on the sums are infinite. This is because a boundary condition 
on the concentration is needed in order to be able to evaluate the square-gradient term at 
the boundary of the container. To this end, it will be assumed that outside the container, 
there is a uniform concentration denoted p^ so that, formally, p i+1 / 2 = Poo f° r * — N. This 
concentration could be zero, in the case of a truly isolated system, or could be set to a 
finite value in order, e.g., to approximate nucleation in an infinite system. As it stands, the 
squared gradient term requires the concentration derivative at intermediate points. This is 
approximated as 
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giving 



OO OO / rj v 2 
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or, with the boundary condition, 

N N+l , _ s 2 

/3F — /3Foo = A ^ 47rrJ +1/2 /3/ (p, +1/2 ) + A £ 2tt^ ( ^ Pj-1/2 
3=0 3=0 ^ ' 

where /Si 7 ^ is the energy of the system with uniform concentration p^. 



3. The stochastic differential equation 

With this discretization, the SDE becomes 



dm i (t) _ n ra»+i (t) -mj-i (t) OF /on ^+l (*) ~ wti-i (*) /I. /.x 1 / ■ ^ A r 

dt ~ 2A 9Am,(t) + V 2A VA €iW ' " " 

(14) 

with (£) £ ■ (t')) = 5 (£ — £') There are only equations for m; for z up to N since, as 
discussed above, the wall is at position N + 1/2 and the total mass in the container, mjv+i, 
is conserved. Using Eq.([5]) and (Q gives 

^ = 4^ 2 a( (1 - M ^T- (1 -^ ) ^J' (15) 

where it the invariance of m^ + i and of mo = have been explicitly indicated. Expanding, 
this gives, 

with 

_1 n r , r i+3/2Pi+3/2+ r i+l/2Pi+l/2 f 1 1 

% = (1 " M ^3 -2 ^ - - 2 S i+lj J (17) 
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and 



Si\r) Si+n — (1 — 5 in) 5i 



9ij = WlA^ r^ 2 Ti+wPi+W + r UliPi-^ 0<iJ<N 

1+1/2 (18) 
It is easy to check that the fluctuation-dissipation relation g" 1 = qq T holds and that the 

discretized equation remains Ito-Stratonovich equivalent^. 

Just as one of the mass variables is fixed by the constraint of constant mass, so one of 
the concentration variables is also redundant. In fact, from the previous equation it follows, 
since the total mass is conserved, that 

dpN+1/2 _ 1 fdm N+1 dm N \ 1 dm N 

(iyj 



dt 47rr 4 2 +1/2 A \ dt dt J 47rr 2 +1/2 A dt 

There are two possibilities: either one proceeds using this as the equation of motion or one 
eliminates the variable Pn + i/2 altogether by means of the expression for the total mass giving 

r N+l/2pN+l/2 = JZ m N+l ~ Yl r i+l/2Pi+l/2 (20) 

i=0 

so that in fact the remaining variables are the masses for 1 < j < N or, equivalently, the 
densities up to /0jv-i/2- When evaluating derivatives, one must then use, e.g., 

d^F d(3F 



dPi+i/2 dPi+1/2 



dpi+i 



/2 



| dpN+i/2 d$F 

dPi+l/2 dPN+l/2 

Pn+i/2 

r 2 i+l/2 df3F 



r 



2 

N+l/ 



1/2 ^PN+I/2 



PN+l/2 

The notation indicates that in the first term on the right hand side the derivatives are 
evaluated holding p^+\/2i as we ^ as a ^ other concentration p J+ i/2 except for j = i, constant. 
Note that the second term on the right plays no role when discretizing the SDE, Eq. ffTB]) . 



4- Stationary solutions 

A special role will be played by the stationary solutions to the deterministic dynamics 
(i.e. the dynamics with no noise). These will include both the metastable initial and final 
states and the metastable critical cluster separating them. They are simply the solutions of 

dF 

= 0, 1 < i < N (22) 

orrii 



under the constraint that the total mass is conserved (so that all derivatives are understood 
to be evaluated with mN + i fixed). Equivalently, in terms of the concentration this gives 

l-StN dF 1 dF 



0, 1< % < N. 



(23) 



47n 7+l/2 9 Pi+l/2 47 ™ll/2 dPi-1/2 

where the derivatives are evaluated as indicated in Eq. (l21l) . so that this could also be written 



as 



dF 



47rrf_ d Pi _ 1/2 
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dF 



^ r N-l/2 d PN-l/2 
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1 dF 



1< % < N- 1 



(24) 
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Finally, by introducing an additional variable, p, these N equations can be replaced by the 
N + 1 equations 



dF 



dp 



t-l/2 



4mj_ 1/2 n, l<i<N 



(25) 



Pj + l/2 



N 



m N+1 = A 4vrr J 2 +1/2 p j+1/2 

3=0 

for the variables p i+1 / 2 ) < % < N and //. In this form, it is clear that the stationary 
solutions simply correspond to extrema of the free energy functional subject to the condition 
of constant mass, as one would expect. A convenient method of solution of these equations 
is given in the Appendix |A] 



C. Nucleation pathway 

In principle, one could simulate the discretized SDE but the goal here is to more di- 
rectly characterize the nucleation pathway. Since the problem of interest is the nucleation 
of a dense phase, with average bulk concentration p h from a dilute phase with average 
bulk concentration p v , one approach to characterizing the nucleation pathway is to look 
for the most likely path^~- connecting an initial state with p(r;t = 0) ~ p v and final state 
p (r;t = T) ~ p x . For an infinite system, these conditions would be strict equivalences but 
for a finite system, there will in general be spatial variation of the concentration, even in 
the initial state due to the boundaries. Furthermore, for a finite container, since mass is 
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conserved, there cannot be a complete conversion from the initial state p (r; t — 0) ~ p v 

to the final concentration, p t : the lowest energy state that conserves mass will consist of a 

droplet inside of which the concentration is roughly that of the final state, p t , and outside 

of which the concentration decays to some value that may or may not be near p v depending 

on the boundary conditions. 

It has previously been shown that in the weak noise limit, in which the thermodynamic 

driving force is large compared to typical noise fluctuations, the most likely path connecting 

the two metastable states will pass through the critical cluster-^. Furthermore, the MLP 

can be determined by following the deterministic equations of motion forward and backward 

in time. In practice, this means that one starts at the critical cluster and then perturbs the 

system slightly, first in the direction of one metastable state and then in the direction of the 

other, and follows the deterministic dynamics given in the present case by 

dm (t) _ n m m (t) -m^i(t) OF 

U 2A dAm, (t) ' iS ^ iV W 

or, equivalently, 

where, again, the derivatives on the right are evaluated using Eq. (121j) . There are several 
important points to be made about this procedure. First, even though there is no noise 
term in Eq.( [26]) . this is a noise driven process: Eq.f l26|) is simply a convenient and exact 
means of determining the most likely path which is obviously a concept rooted in the fact 
that the full dynamics includes fluctuations^. Second, the fact that one follows the path in 
both directions away from the critical cluster, or equivalently, that one follows the dynamics 
forward and backwards in time, does not in any way imply time-reversal invariance. In 
fact, the underlying dynamics is dissipative and so is definitively not time- reversal invariant. 
Again, the fact that this procedure serves to determine the MLP is a mathematical result 
that does not imply anything concerning the time-reversal properties of the actual stochastic 
dynamics. Finally, the first part of the prescription, namely perturbing away from the critical 
cluster, must be made more precise. One must perturb away from the initial state because 
the deterministic driving force in the critical cluster is zero: it is an unstable stationary point 
of the dynamics. Thus, in principle one must move an infinitesimal distance away from it 
to initiate the dynamics, but in practise one must of course make a finite perturbation with 
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the understanding that the final results are approximate and approach the exact result only 
as the size of the initial perturbation goes to zero. To construct the perturbation, one first 
solves a generalized eigenvalue problem as described in Appendix [B] which typically results 
in the identification of a single, unstable eigenvalue. The associated eigenvector, Vi, defines 
the unstable direction and the initial state is simply 

p l+1/2 = p* +1/2 ±ev u 0<i<N-l (28) 

where e is a parameter controlling the size of the imposed perturbation. 

The next step is the numerical integration of Eq. (l27|) . Because of the Ito-Stratonovich 
equivalence of the (spatially-) discretized SDE there is no constraint on the discretization 
scheme used for the time variable. For example, with discrete time variables = it, the left 
hand side is simply (^4+1/2(^+1) — Pi+i/2^j))l T while the right hand side can be evaluated 
at tj, giving an Euler scheme, or at tj + r/2 giving a more efficient implicit scheme. Further 
details are given in Appendix 

Finally, it is convenient to introduce a parameterization for the nucleation pathway. As 
described elsewhere^, the structure of the stochastic model induces a natural measure of 
distance in concentration space between the concentration at time t\ and that at time t 2 
given by 

By definition, this distance increases monotonically along any proposed nucleation pathway 
and therefore serves as a natural, induced reaction coordinate. 



III. EXAMPLE: NUCLEATION OF DENSE PROTEIN SOLUTION FROM 
DILUTE SOLUTION 

In order to illustrate the theory developed above, detailed calculations for a model glob- 
ular protein^ have been performed. The effect of the solvent was approximated, crudely, 
by assuming Brownian dynamics of the (large) solute molecules which also experience an 
effective pair interaction for which the ten Wolde-Frenkel interaction potential 

00, r < a 

v(r)={ 4e ( f 1 X 6 ( , \*\ . (30) 
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FIG. 1. The concentration distribution in the metastable, nearly uniform initial state (solid line), 
the critical cluster (broken line) and the final stable cluster (dash-dot line) as a function of distance 
from the center of the cluster. 

was used with a = 50 which is then cutoff at r c = 2.5a and shifted so that v (r c ) = 0. The 
temperature is fixed at fc_gT = 0.375e and the equation of state computed using thermody- 
namic perturbation theory. The transition studied was that between the dilute phase and 
the dense protein phase which, in the present simplified picture, is completely analogous 
to the vapor-liquid transition for particles interacting under the given pair potential. The 
gradient coefficient, K, is calculated from the pair potential using the approximation given 



in Ref. 
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f3K~ d 5 Pv(d) + — (2d 2 - 5r 2 ) f3v (r) r 2 dr (31) 

45 15 J d 



where d is the effective hard-sphere diameter for which the Barker-Henderson approximation 
was usecr^. For the temperature used here it was found that f$K = 1.80322cr 5 . 

The nucleation pathway has been evaluated for a spherical box with a radius or 20.63cr 
using a lattice of 160 points with a spacing equal to 0.126a and with the external concentra- 
tion set to zero. The time step used was lO -4 /}" 1 ^ 2 . The only remaining control parameter 
is the initial mass of material in the container. For the chosen temperature, the coexistence 
densities of the dilute and dense phases are pa 3 = 0.071 and per 3 = 0.663, respectively. The 
initial condition for the calculations presented here was a uniform concentration which was 
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1.5 times the dilute phase coexistence concentration which serves to fix the total mass. The 
various metastable states were then determined under the constraint that the total mass 
does not vary. Figure [T] shows the initial metastable state, the critical cluster and the final, 
stable state. Since the material in the center of the cluster that eventually forms is near the 
bulk state, it is clear that the if the initial state contains too little matter, no cluster will 
form. In fact, if the amount of matter in the system is scaled with the radius, then as the 
radius of the system is reduced, the final stable cluster has smaller and smaller radius until 
at some point it is smaller than the critical cluster: i.e., there is no barrier separating the 
initial and final states. 

One quantity of particular interest is the size of the cluster. In previous work on infinite 
systems, this was characterized by both the equimolar radius and by identifying the point 
at which the difference of the concentration from the background was sufficiently small. In 
the present case of relatively small, finite systems that conserve mass, the equimolar radius 
is not very useful so only the second measure is employed. Even then, the "background" 
concentration has to be reinterpreted: here, it is taken to be the concentration a distance of 
la from the wall (so as to avoid boundary effects), a point denoted as r& since it defines the 
"background". The radius is then determined by starting at the center of the cluster, i.e. 
at r = 0, and increasing r until a point R is found at which 

p(R)=p(r b ) + e(p(0)-p(r b )) 

for some prescribed value of e, typically either e = | (which will be referred to as the half 
width) or e = 0.1 (which will be referred to as the total width). 

The determination of the two branches of the pathway was initiated by perturbing in the 
direction of the unstable eigenvector with the proportionality constant being adjusted so as 
to produce a change in free energy of AF = 0.005ksT from the critical cluster. Figure [2] 
shows the evolution of the cluster radius, the central concentration and the free energy differ- 
ence as functions of distance as well as of the cluster size along the nucleation pathway. The 
results are similar to those obtained previously for an infinite system using parameterized 
concentration profiles. In particular, the fact that the process of nucleation involves two 
separate steps: first, a long wavelength, small amplitude concentration fluctuation forms. 
This stage is characterized by a large initial radius that decreases while the central concen- 
tration increases. At a certain moment, when the concentration is near its bulk value, the 
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FIG. 2. The left panel shows the free energy difference from the initial state, AF/AF C , the con- 
centration at the center of the cluster, p(0)/p c (0), and the cluster radius, Rhai f I (4c) , as functions 
of distance along the nucleation pathway. The energy and concentration are scaled by their values 
in the critical cluster whereas the radius has been arbitrarily scaled to 4<r. The vertical line shows 
the point at which the radius reaches its minimum. The right panel shows two measures of the 
cluster size as a function of distance along the pathway. 

radius begins to increase and, in this second stage, the evolution of the cluster follows the 
classical scenario of monotonic growth of radius and mass at constant, near-bulk interior 
concentration. 

There is one significant difference in the pre- and post-critical cluster growth as shown in 
Fig. [3j Growth of the post-critical cluster is deterministic and driven by the fact that the free 
energy of the system decreases with increasing cluster size so that material is actively drawn 
into the cluster. The fact that the low-concentration protein drawn into the cluster can only 
be replaced by the diffusion of new material from further away results in the formation of 
a region with a lower concentration of protein than in the mother phase, further from the 
cluster; in other words, in the formation of a depletion zone. The formation of a depletion 
zone is expected in classical treatments of cluster growth^. The process is governed by 
the fact that the lower the concentration in the mother phase, the larger the corresponding 
critical cluster. If material is drawn into the growing cluster faster than new material can 
replace it by diffusion, (a situation that is always assumed to occur for sufficiently large 
clusters since the thermodynamic driving force increases with cluster size), then the cluster 
will lower the concentration in the adjacent region to the point that the thermodynamic 
driving force is sufficiently small that the rate of addition of material to the cluster will be 
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FIG. 3. The evolution of the concentration along the nucleation pathway. The left panel shows 
the early stages of formation of the pre-critical cluster (increasing concentration at r = corre- 
sponds to increasing distance along the nucleation pathway). The large spatial extent of the initial 
concentration perturbation is evident and the inset shows the monotonic variation of the concen- 
tration as a function of distance from the center. The right panel shows part of the concentration 
distribution during post-critical growth (the maximum concentration moves to increasing radii as 
one progresses along the nucleation pathway). The critical cluster, left-most curve, and final state, 
right-most curve, show no depletion zone, but the presence of a depletion zone in intermediate 
clusters is clear. 

equal to the rate of replenishment by diffusion. Assuming diffusion is slow, this will force 
the cluster to be very near the critical cluster for the adjacent depleted mother phase. 

In contrast to the formation of the depletion zone for post-critical clusters, Fig. [3] shows 
that no depletion zone exists for pre-critical clusters. The reasoning is the opposite to that 
for post-critical growth: the pre-critical cluster is unstable and its growth is governed by 
unlikely fluctuations. If the concentration outside the cluster is unnecessarily low, then the 
cluster will be even more unstable (since, as pointed out above, lower concentration in the 
mother phase increases the size of the critical cluster). Hence, the most favorable way to 
form an unstable cluster is one that does not create any region of lowered concentration. 
Of course, since mass is conserved, the concentration outside the cluster must decrease as 
the cluster grows, but, as can be seen in the figure, the concentration is at a maximum just 
"outside" the cluster (however that is defined) and decreases monotonically as one moves 
further away from it so that the concentration of the mother phase just outside the cluster is 
as large as possible, i.e. the maximum of the concentration anywhere in the mother phase. 
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This illustrates a point made previously^, namely that the initial formation of a cluster as 
a localized object with small radius would necessarily involve the formation of a depletion 
zone, due to mass conservation, and this would clearly be unfavorable since it would make 
any (unstable) pre-critical cluster even more unstable. Instead, what the present results 
show is that the most likely nucleation pathway involves a long-wavelength concentration 
fluctuation that contains and localizes enough mass to form the critical cluster. 

IV. CONCLUSIONS 

A discretized version of the stochastic field equations for fluctuations in a fluid of Brown- 
ian particles in the over-damped limit have been given. The most likely path for nucleation 
was computed in the weak-noise limit and it was found that the pathway obeyed the same 
two-stage mechanism as previously found for infinite systems. Namely, nucleation begins 
with a long-wavelength, small- amplitude fluctuation that evolves by becoming spatially more 
localized until the concentration inside the fluctuation reaches near-bulk values. This is fol- 
lowed by a second stage of spatial growth of the cluster that is similar to that assumed in 
CNT. 

It was observed that an important difference between pre-critical and post-critical growth 
is that in the former case, no depletion zone forms whereas in the latter, a depletion zone 
is clearly in evidence. The existence of depletion zone is not surprising and, indeed, is 
assumed in classical treatments of cluster growth. However, its presence is important in 
distinguishing the results of the present, dynamical approach to nucleation from those of 
other approaches. One alternative method is the use of equilibrium DFT. The determination 
of the properties of the critical cluster via DFT is relatively uncontroversial, however the 
use of DFT to determine the pathway, while widely practiced, is harder to justify. This 
is because DFT is a strictly equilibrium theory^ whereas nucleation is, by its nature, a 
nonequilibrium, fluctuation-driven processes^. It is therefore unclear why DFT should be 
the right tool to describe the nucleation pathway. This is reflected in the fact that it is 
difficult to map out the pathway using DFT which is based on minimizing a free energy 
functional. The problem is that non-critical clusters are obviously not stationary points of 
the free energy so that one can only determine their properties using DFT if some sort of 
external constraint is applied^— so as to stabilize them. This leads to the problem that 
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there is no unique constraint and different constraints can lead to different results^ 1 ^. 

In answer to this problem, another approach has gained popularity in recent years. It 
involves the determination of the minimum free energy pathway connecting the critical 
cluster to the initial and final states^ - — . While this does not involve the imposition of 
any constraints, it nevertheless suffers from a similar arbitrariness as the MFEP can only 
be calculated once a method for determining distances in concentration-space has been 
adopted^. Again, without further input the choice of metric in concentration space is 
arbitrary and the method ends up being no better founded than the older constraint-based 
methods. 

Aside from these criticisms of the internal logic of the DFT methods, the present results 
show a further qualitative failing that demonstrates that free energy considerations, without 
dynamics, cannot be enough. As discussed above, the formation of a depletion zone is an 
inevitably consequence of the fact that the thermodynamic driving force for the incorporation 
of new material into the cluster grows as the cluster grows and must, at some point, cause 
material to be drawn into the cluster faster than it can be diffusively replaced. None of 
the DFT studies report seeing a depletion zone for super critical clusters. In all cases, the 
concentration seems to decay monotonically to the surrounding bulk. The is a clear sign 
that approaches that do not include dynamics are not sufficient. 

An alternative to DFT is its dynamical cousin, Dynamical DFT (DDFT)^£ The equa- 
tions for DDFT are in fact the same as those given above, Eq.(TjQ), if the noise term is 
dropped or Eq.Q if it is understood that only forward evolution in time is relevant (for 
DDFT). The latter point is important: DDFT can only describe a system falling down a 
free energy surface: by its nature, since it involves no fluctuations, DDFT cannot describe 
barrier crossing. To overcome this, some studies^ have relied on the inclusion of fluctuations 
in the initial condition for a DDFT calculation but at best this describes a physical system 
perturbed by an initial, external shock. The results will necessarily be dependent on, e.g., 
the amplitude and spectrum of the noise and thus is subject to a similar arbitrarity to that 
of static DFT and for the same reason: both theories are based on equilibrium physics and 
so cannot naturally describe nonequilibrium states. On the other hand, DDFT does have 
the virtue that the post-critical growth is correctly described and it will certainly correctly 
describe the depletion zone. 

The dynamical theory presented here and elsewhere does, in some sense, provide an ex- 

1(3 



tension of DDFT to barrier crossing problems since the MLP is described by DDFT-like 
equations. However, there are still conceptual differences that separate the two approaches. 
The most important is the nature of the "free energy" that enters into the dynamical equa- 
tions. In DDFT, the free energy typically enters via a local equilibrium approximation. 
Because of the context (i.e. that one is dealing with ensemble-averaged quantities) this is 
the true, macroscopic free energy functional. In the dynamical approach, the status of the 
"free energy" is more problematic. It is certainly not the macroscopic free energy. In fact, 
let us write the DDFT equations in the form 

\ ^ / p(r)^{p(r ;i )> 

where the brackets (...) indicate an ensemble average and where F [p] is the equilibrium DFT 
functional. On the other hand, the noise-averaged stochastic equations, Eq.([T]), are 

9M»_ BV // pM)v ^|Pl 



» \\ VMr)/ p(rH( , (r;t)l 

where the double brackets ((•••)) indicate an average over the noise. Technically, one cannot 
identify the ensemble- averaged concentration, (p(r;t)), with the noise-averaged concentra- 
tion ((p(r;t))) because the noise-averaged concentration also involves a coarse-graining. 
However, if the spatial variations are sufficiently weak, then it is reasonable to expect these 
quantities will be more or less the same and so one expects the right hand sides of these 
equations to become the same. It is then seen that equality of the right hand sides depends 
on a series of factorization approximations, 



5F[p] 



5p(r) 



<(pM)»v 

p(rH«p(r;t))) 

Interestingly, the assumption of weak noise, which underlies the development of the dynam- 
ical approach, might be expected to justify this series of approximations so that a certain 
consistency between the dynamical theory and DDFT might be seen to emerge. In any case, 
it is worth noting that in the theory of critical phenomena, it is generally accepted that a 
mean field approximation of the form hard-core plus mean-field treatment of the attractive 
part of the potential is a reasonable approximation for the coarse-grained free energy^. 
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Since this is at present the dominant model used in DFT calculations^, it might be argued 
that the use of this model for the coarse-grained free energy is actually more justifiable than 
is its use in DFT calculations. 
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Appendix A: Solution for the critical cluster 

With the discretized free energy, Eq. (fl3]) . the equations for the critical cluster, Eq.( l25|) . 
become 

2 ( P3/2 ~ Pl/2 



r\ n M' (ft /2 ) + Kr\ J = - V < A1 > 



r N+l/2p 



and 

N 



m N+l = 4tt r* +1/2 p i+1/2 . (A2) 

i=0 

To solve these, one begins with a guess, say p i+1 / 2 , and denote the actual solution by p* +1 / 2 - 
Defining 

x i = Pi+1/2 - Pi+1/2, <i < N (A3) 

XN+l = fi* — fi 

and 

Vi = Pf (Pi+1/2) + 2 ^ A ( Pi +V2 - ft-l/2) - r i+l (pi+3/2 - ft+l/2)) > < i < N - 1 

r i+l/2^ 



(A4) 



VN = (3f (PiV+1/2) + "2 (PjV+l/2 - PN-l/2) ~ r N+l (Poo - PN+l/2)) 



r N+l/2 



N 



y N+1 = 47T Y r i+\/2Pt+\/2 ~ m N+l 



i=0 
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and expanding to first order gives 

^ V^+l/2/ \ r i+l/2 A / 

K ( r- i \ 2 

- (1 - 5 iA r) — — x m - y N+ i = p- yi, < % < N (A5) 

A V r i+l/2/ 

and 

4tt r 2 l+1/2 x i+1/2 = y N +i- (A6) 

i=0 

These equations are tridiagonal except for the last row and column and can be efficiently 
solved by a trivial modification of the Thomas algorithm (see, e.g., Ref.— ). Thus, the solution 
to these equations is found by starting with the boundary conditions, p^ and m^ +1 , and an 
initial guess, p i+ i/ 2 and p (for the latter, a reasonable guess is /' (p^)) and then iterating 
this linear system to stability. Alternatively, one can simply choose a value for p, eliminate 
Un+i as a variable (i.e. set i/n+i = in Eq. flA5j) ) and drop Eq. (1A6|) and solve for the values 
of Xj. The value of the mass is then evaluated post hoc. 



Appendix B: The unstable direction 

To determine the unstable direction, we first note that the deterministic dynamics 
given in Eq. (f27|) can be interpreted as being gradient- driven motion in a non- Euclidean 
geometry^^ 1 ^. As such, the unstable direction is determined from the generalized eigen- 
value problem 

N-l ^ N-l 

= A$>^, 0<.<A-1 (Bl) 

j=0 d Pi+l/2°Pj+l/2 j=Q 

and for convenience we write 



with 



and 



Rij = 5ijr i+1/2 (B3) 

9i/ = (1 - M (ri +3 /2Pi+3/2 + r i+l/2Pi+l/2) ~ $i+lj) ( B4 ) 

- (1 - 6 m ) (r- +1/2 A+i/2 + r Li/2Pi-i/2) (h-ij - > < i, j < N - 1 
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which we also write as 



9ij = - (diSji-i + ajdji+i) + ((1 - 6 i0 ) di + a i+ i) <% < i, j < N - 1 (B5) 

where 



a; = r, 



" +l/2Pt+l/2 + r i~l/2Pi-l/2, <ij < N -1 (B6) 



AT-l 

&i = 



The inverse of g 1 is easily found to be 

E — ■ ( B7 ) 

l=max(j,j) 

Thus completing specification of the generalized eigenvalue problem. 

Appendix C: Integration scheme 

The deterministic equations can be written as 

= Vi+1/2 (Pi-3/27 Pi-1/2) Pi+1/27 Pi+3/27 Pi+5/2) (CI) 



with 



2P+1/2 (^Pi-3/2 5 Pi-1/2 5 Pi+1/2 5 Pi+3/2 5 Pi+5/2) 

(C2) 

v ^+3/2^+3/2 + ^+l/ 2 A+l/2 / 1 <9/3F 1 <9/3F 

87rr '+l/2 A3 ^+1/2^+1/2 ^+3/2^+3/2^ 

T.- 1 1 in Da l 1 /o H 

+ (1 - M 



r ^i/2Pi+i/2 + r?_ 1/2 ft_i /2 ( 1 dfiF 1 



^ r '+l/2 A3 V »i_i /2 ^-V2 r '+l/2 <^ 



i+1/2 



and 



= A4?r f Z 3 /' fe+1/2) + r 2 K A 2 [- r 'P,-i/2 + (*i + r^i) P m/2 - r? +1 p i+3/2 ] 

r i+l/2 °Pi+l/2 \ r i+l/2 lA 

(C3) 

Discretizing in time with time step 8 t gives a simple Euler scheme, 

- P t\/2 = fclffi/a (G4) 

with 

(») f» » » » » 



2/i+l/2 - 2/*+l/2 ^Pj-3/2' Pi-1/2' Pi+1/2' Pi+3/2' Pj+5/2 J 
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Unfortunately, this simple scheme requires a very small time step and is therefore inefficient. 
As is often the case^, the integration can be made much more efficient by introducing an 
implicit scheme. Define 

„("+!) _ >+!) _ >) 
and evaluate the source term at the temporal "midpoint" using 

(n+1/2) = 1 f (n+1) („) \ _ (n+1) _ 1 (n+1) rm 
^i+1/2 — 2 V^+l/2 ^ Pi+1/2) ~ Pi+l/2 2 1+1 1 2 ^ ' 

(n+1) 



and expanding to first order in x^ n+ ' to get 



2 

(n+1/2) _ (n) T (n) (n+1) 

y?+i/2 — yj+i/2 + »i j+i/2 l^'J 
i=-2 



with Tjj = Tjj and 



8vrA3 ^ +1/2 «9p m/2 rf +3/2 9p i+3/2 / 
1 / 1 d/3F 1 <9/3F 



+ (1 - S i0 ) 



87rA 3 V»*i-l/2 9 Pi-l/2 ^+1/2^+1/2 



2 i n ) j_ 2 ( n ) / 

fl A \ ^+3/2Pi+3/2 + r i+l/2Pi+l/2 f 1 R fll ( in) i 

- (i - ^/ ^ +1/2 J H 

2 ( n ) I 2 ( ra ) / 

, , , r i+l/2Pi+l/2 + r i-l/2P»-l/2 [ 1 nr/ffM 

-(l-*o) 2^2^ l^fA/ (^ +1/2 




and 



87rr i+l/2^ \ r i+l/2^Pi+l/2 r i+3/2^Pi+3/2 / 

2 ( n ) i 2 ( n ) / / 9 

r i+3/2/^+3/2 r i+l/2Pj+l/2 / 1 „ / ( n ) 

i+2 



+ {1 _ Sn) m ^ ^ ^ + _ ^ + ^ 

2 ( n ) I 2 ( n ) / 2 

. /-, c \ r i+l/2Pi+l/2 + r i-l/2Pi-l/2 A / r i+1 

+ CI - M ^ p 

Z 'j+l/2 ^ \'i+l/2 > 

and 

„2 



rr +3/2Pi+3/2 + r i+l/2Pi+l/2 K ( r i+2 

r 

The integration scheme then becomes 



( n ) 2 (n) 

! + r i+l/2Pj+l 

2 

1+1/2 V'i+3/2 



!+^ = 6tv$! 1/2 + ** £ ^ (cii) 

i=-2 
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Because of the banded nature of these equations, they can be solved quite efficiently and 
the resulting integration scheme is much more stable than the simple Euler scheme. 
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